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1 The eigenvalue distribution function 

For a.n. N X N matrix An, the eigenvalue distribution function ^ (c.d.f.) F^" {x) is defined as 

pA^ ^ ^ Number of eigenvalues of Ajy < x 

As defined, the e.d.f. is right continuous and possibly atomic i.e. with step discontinuities at discrete 
points. In practical terms, the derivative of (1), referred to as the (eigenvalue) level density, is simply the 
appropriately normalized histogram of the eigenvalues of Af^ . The MATLAB code histn we distributed earlier 
approximates this density. 

A surprising result in infinite RMT is that for some matrix ensembles, the expectation i?[i^^"(a;)] has a 
well defined i.e. not zero and not infinite limit. We drop the notational dependence on in (1) by defining 
the limiting e.d.f. as 

F^{x) = lim E{F^'^[x)\. (2) 

This limiting e.d.f. ^ is also sometimes referred to in literature as the integrated density of states [2, 3]. Its 
derivative is referred to as the level density in physics literature [4]. The region of support associated with 
this limiting density is simply the region where dF^{x) ^ 0. When discussing the limiting e.d.f. we shall 
often distinguish between, its atomic and non-atomic components. 



2 The Stieltjes transform representation 

One step removed from the e.d.f. is the Stieltjes transform which has proved to be an efficient tool for 
determining this limiting density. For all non-real z the Stieltjes (or Cauchy) transform of the probability 
measure F-^{x) is given by 

mA^z) = ( -^dF^{x) Im z^O. (3) 



The integral above is over the whole or some subset of the real axis since for the matrices of interest, 
such as the Hermitian or real symmetric matrices, the eigenvalues are real. When we refer to the "Stieltjes 
transform of A" in this paper, we are referring to to^(z) defined as in (3) expressed in terms of the limiting 
density dF^[x) of the random matrix ensemble A. 



^This is also referred to in literature as the empirical distribution function [1]. 

^Unless we state otherwise any reference to an e.d.f. or the level density, in this paper will refer to the corresponding limiting 
e.d.f. or density respectively. 

^While the Stieltjes integral is over the positive real axis, the Cauchy integral is more general [5] and can include complex 
contours as well. This distinction is irrelevant for several practical classes of matrices, such as the sample covariance matrices, 
where all of the eigenvalues are non-negative. Nonetheless, throughout this paper, (3) will be referred to as the Stieltjes 
transform with the implicit assumption that the integral is over the entire real axis. 



The Stieltjes transform in (3) can also be interpreted as an expectation with respect to the measure 
F^{x) such that 



mAiz) = Ex . (4) 

\_x — z 

Since there is a one-to-one correspondence between the probabihty measure F^{x) and the Stiehjes 
transform, convergence of the Stieltjes transform can be used to show the convergence of the probability 
measure F^{x). Once this convergence has been established, the Stieltjes transform can be used to yield 
the density using the so-called Stieltjes-Perron inversion formula [6] 

dF^{x) 1 . _^ 

= — lim Im mA[x + ic,)- (5) 

dx TT ►o 

When studying the limiting distribution of large random matrices, the Stieltjes transform has proved to 
be particularly relevant because of its correspondence with the matrix resolvent. The trace of the matrix 
resolvent, Ma{z), defined as Ma{z) = {Aj^ — zl)~^ can be written as 



^ 1 

tr[M^(.)]^^^ (6) 

i—l 

where A; for i = 1,2, . . . , N are the eigenvalues of A]\j. For any An, Ma{z) is a non-random quantity. 
However, when yl^r is a large random matrix, 

m^(z)= lim l-tr[MA{z)]. (7) 

The Stieltjes transform and its resolvent form in (7) are intimately linked to the classical moment problem 
[6]. This connection can be observed by noting that the integral in (3) can be expressed as an analytic 
"multipolc" series expansion about z = oo such that 

/oo ^. oo p f, 

Y.^dF^i-)-~j:^dF^i-) 

k=0 k=0 /„N 
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where = J x'^dF^{x) is the fc*'' moment of x on the probability measure dF^{x). The analyticity of 
the Stieltjes transform about z = oo expressed in (8) is a consequence of our implicit assumption that the 
region of support for the limiting density dF^{x) is bounded i.e. limj^^oo dF^{x) = 0. 

Incidentally, the 77-transform introduced by Tulino and Verdii in [7] can be expressed in terms of m{z) 
and permits a scries expansion about z = 0. 

Given the relationship in (7), it is worth noting that the matrix moment is simply 



= lim ^iv[A%] = / x'^dF' 

N—>oo iV J 



ix). (9) 



Equation (8), written as a multipole series, suggests that a way of computing the density would be to 
determine the moments of the random matrix as in (9), and then invert the Stieltjes transform using (5). For 
the famous semi-circular law, Wigncr actually used a moment based approach in [8] to determine the density 
for the standard Wigncr matrix though he did not explicitly invert the Stieltjes transform as we suggested 
above. As the reader may imagine, such a moment based approach is not particularly useful for more general 
classes of random matrices. We discuss a more relevant and practically useful Stieltjes transform based 
approach next. 



3 Stieltjes transform based approach 



Instead of obtaining the Stieltjes transform directly, the so-called Stieltjes transform approach relies instead 
on finding a canonical equation that the Stieltjes transform satisfies. The Marcenko-Pastur theorem [9] was 
the first and remains most famous example of such an approach. We include a statement of its theorem in 
the form found in literature. We encourage you to write this theorem, as an exericse, in a simpler manner. 

Theorem 1 (The Marcenko-Pastur Theorem ). Consider an N x N matrix, Bm- Assume that 

1. Xn is an nx N matrix such that the matrix elements X^^ are independent identically distributed (i.i.d.) 
complex random variables with mean zero and variance 1 i.e. X"j £ C, E[X^j] ~ and = 1. 

2. n = n{N) with n/N c> as N ^ oo. 

3. Tn ~ diag{T" , T2 , . . . , t") where r" e K, and the e.d.f. 0/ {t", . . . , r"} converges almost surely in 
distribution to a probability distribution function H[t) as N 00. 

4. Bn = An + jjX*TnXn, where An is a Hermitian N x N matrix for which F^'^ converges vaguely to A 
almost surely, A being a possibly defective (i.e. with discontinuities) nonrandom distribution function. 

5. Xn, Tn, and Ajsi are independent. 

Then, almost surely, F^'^ converges vaguely, almost surely, as N ^ 00 to a nonrandom d.f. F^ whose 
Stieltjes transform m{z), z £ C, satisfies the canonical equation 

We now illustrate the use of this theorem with a representative example. This example will help us 
highlight issues that will be of pedagogical interest throughout this semester. 

Suppose An ~ i.e. Bn = j^X*TnXn. The Stieltjes transform of An, by the definition in (3), is then 
simply 

m^(z) = -i- = -i. (11) 

— z z 

Hence, using the Marcenko-Pastur theorem as expressed in (10), the Stieltjes transform m(z) of Bn is given 

by 

"^("^^--—— TTZhW (12) 

J l+rm(^) 

Rearranging the terms in this equation and using m instead of m{z) for notational convenience, we get 

'•^f^'^m. (13) 



TO J 1 + rm 

Equation (13) expresses the dependence between the Stieltjes transform variable to and probability space 
variable z. Such a dependence, expressed explicitly in terms of dHir), will be referred to throughout this 
paper as a canonical equation. Equation (13) can also be interpreted as the expression for the functional 
inverse of to(z). 

To determine the density of Bn by using the inversion formula in (5) we need to first solve (13) for m(z). 
In order to obtain an equation in m and z we need to first know dH[T) in (13). In theory, dHir) could 
be any density that satisfies the conditions of the Marcenko-Pastur theorem. However, as we shall shortly 
recognize, for an arbitrary distribution, it might not be possible to obtain an analytical or even an easy 
numerical solution for the density On the other hand, for some specific distributions of dHij), it will indeed 
be possible to analytically obtain the density Wc consider one such distribution below. 

Suppose Tn = I i.e. the diagonal elements of r„ are non-random with d.f. dH[T) = 5[t — 1). Equation 
(13) then becomes 

z^-- + -^. (14) 
m 1 + m 



Rearranging the terms in the above equation we get 



zm{l + m) = —(1 + m) + cm (15) 

which, with a bit of algebra, can be written as 

TO^z + m(l - c+ 2) + 1 = 0. (16) 

Equation (16) is a polynomial equation in m whose coefficients arc polynomials in z. We will refer to such 
polynomials, often derived from canonical equations as Stieltjes (transform) polynomials for the remainder 
of this paper. 

As discussed, to obtain the density using (5) we need to first solve (16) for m in terms of z. Since, from 
(16), we have a second degree polynomial in m it is indeed possible to analytically solve for its roots and 
obtain the density. 




Figure 1: Level density for Bn = {l/N)X*Xn with c = 2. 
This level density, sometimes referred to in the literature as the Marcenko-Pastur distribution, is given 



by 



^— = max 0, 1 - cj<5(.) + — /[ ,_,,^] (17) 



where b± = (l±-y/c)^ and I[b_^b+\ is the indicator function that is equal to 1 for b- < z < b+ and elsewhere. 

Figure 1 compares the histogram of the eigenvalues of 1000 realizations of the matrix = j^X^Xn with 
= 100 and n = 200 with the solid line indicating the theoretical density given by (17) for c = n/N = 2. 
We now consider a modification to the Marcenko-Pastur theorem that is motivated by the sample covari- 
ance matrices that appear often in array processing applications. 



3.1 The Sample Covariance Matrix 

In the previous section we used the Marcenko-Pastur theorem to examine the density of a class of random 
matrices Bn = ■^Ar*T„X„. Suppose we defined the N x n matrix, Yn ~ X*^Tn^'^ , then Bn may be written 
as 

Bn = ^YnY;j. (18) 

1/2 

Recall that T„ was assumed to be diagonal and non- negative definite so r„ can be constructed uniquely up 
to the sign. If Y/v were to be interpreted as a matrix of observations, then Bn written as (18) is reminiscent 



c, 



Figure 2: The matrices i?„ and Cn when n > N. 



of sample covariance matrices that appear in many engineering and statistical applications. However, among 
other things, it is subtly different because of the normalization of the YatY^ by the number of rows N of Y/v 
rather than by the number of columns n. Hence, we need to come up with a definition of a sample covariance 
matrix that mirrors the manner in which it is used in practical applications. 

With engineering, particularly signal processing applications in mind, wc introduce the n x N matrix, 
On = Tn'^Xn and define 

Cn - ^0„0: (19) 

to be the sample covariance matrix (SCM). By comparing (18) and (19) while recalling the definitions of Y^r 
and 0„, it is clear that the eigenvalues of i?„ are related to the eigenvalues of C„. For Bm of the form in (18), 
the Marccnko-Pastur theorem can be used to obtain the canonical equation for msiz) given by (13). Recall, 
that by nisiz) wc mean the Sticltjcs transform associated with the limiting c.d.f. F^{x) of Bn as N —> oo. 
There is however, an exact relationship between the non-limiting e.d.f.'s F^'^{x) and the F'-^"'{x) and hence 
the corresponding Sticltjcs transforms mB~^{z) and mc„{z) respectively. Wc exploit this relationship below 
to derive the canonical equation for C„ from the canonical equation for Bn given in (13). 

Figure 2 schematically depicts C„ and Bn when n > N i.e. when c > 1. In this case, C„, as denoted 
in the figure, will have n — N zero eigenvalues. The other TV eigenvalues of C„ will, however, be identically 
equal to the eigenvalues of Bat. Hence, the e.d.f. of C„ can be exactly expressed in terms of the e.d.f. of 
Bn as 

F^'H-)={^^)ho.oo] + ^F^^{-) (20) 

= (c-l)/(o,oo]+c^^''"(a:). (21) 

Recalling the definition of the Stieltjes transform in (3), this implies that the Stieltjes transform mc„iz) of 
Cn is related to the Stieltjes transform mB^{z) of Bn by the expression 

c - 1 

"^c„ {z) = h CTOBjv {z). (22) 

z 

Similarly, Figure 3 schematically depicts C„ and Bm when n < N i.e. c < 1. In this case, Bn, as denoted 
in the figure, will have N — n zero eigenvalues. The other n eigenvalues of Bn will, however, be identically 
equal to the n eigenvalues of C„. Hence, as before, the e.d.f. of B^ can be exactly expressed in terms of the 
e.d.f. of Cn as 

\ n / n 
^-l)Wl + ^^''"(^) (24) 

Once again, recalling the definition of the Stieltjes transform in (3), this implies that the Stieltjes transform 
mc„{z) of Cn is related to the Stieltjes transform ms^iz) of Bn by the expression 

mB^{z) = -(--l]- + -mc„iz)- (25) 
V c J z c 
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Figure 3: The matrices i?„ and Cn when n < N. 



Muhiplying both sides of (25) by c, we get 



cmB^iz) = -(1 - c)- + mc„{z) 
z 



(26) 



which upon rearranging the terms is precisely (22). Equation (22) is an exact expression relating the 
Stieltjes transforms of C„ and Bjy for all n and A^. As iV ^ oo, the Marcenko-Pastur theorem states that 
tob„(z) niBiz) which implies that the limiting Stieltjes transform mc(z), for Cn can be written in terms 
of mB[z) using (22) as 

c - 1 

mc{z) = hcmB(z). (27) 

For our purpose of getting the canonical equation for mc{z) using the canonical equation for msiz), it 
is more useful to rearrange the terms in (27) and express mB{z) can be written in terms of mc{z). This 
relationship is simply 

mB{z) = -(--l]- + -mciz). (28) 
\c ) z c 

Hence, to obtain the canonical equation for mc{z) we simply have to substitute the expression for mB{z) 
in (28) into (13). With some fairly straightforward algebra, that we shall omit here, it can be verified that 
mc{z) is the solution to the canonical equation 



mc{z) 



dH{T) 



{(1 — c — c z m.c{zy\T — z 



(29) 



nl/2 



Incidentally, (29) was first derived by Silverstein in [10]. He noted that the eigenvalues of -^tJ XnX*T, 
were the same as those of ■^X„X*T„ so that (29) was the canonical equation for this class of matrices as 
well. Additionally, in their proof of the Marccnko-Pastur theorem in [11], Bai and Silverstein dropped the 
restriction on r„ being diagonal so that T„ could be any matrix whose e.d.f. F^'^ H . As the reader may 
appreciate, this broadens the class of matrices for which this theorem may be applied. 



3.2 The (generalized) Wishart matrix 

Revisiting the previous example, when T„ = /, C„ ~ ^XnX*. This matrix C„ is the generalized version of 
the famous Wishart matrix ensemble first studied by Wishart in 1928 [12]. In physics literature, the Wishart 
matrix is also referred to as the Laguerre Ensemble [13]. Strictly speaking, C„ is referred to as a Wishart 
matrix only when the elements of Xn are i.i.d. Gaussian random variables. The canonical equation for C„ 
in (29) becomes 

m = i ^ (30) 

(1 — c — c z m) — z 

which upon rearranging yields the Stieltjes polynomial 



czm^ — (1 — c — z)m +1=0. 



(31) 



Figure 4: The density of the (generahzed) Wishart matrix {M^)i(c)} for different vahies of c. 



Once again, (31) is a second degree polynomial in m whose coefficients are polynomials in z. Hence, as 
before, (31) can be solved analytically to yield a solution for m in terms of z. The inversion formula in (5) 
can then be used to obtain the limiting density for C„. This is simply 



M- = max I 0, 1 - - U(a;) + ^ U[ (32) 



dx \ c J 2ttxc 

where, as before, h± = (l±-y/c)^. As (32) suggests, the limiting density for C„ depends only on the parameter 
c. Hence, for the remainder of this paper, we will use the notation {Ty(c)} to denote the (generalized) Wishart 
matrix ensemble defined as W{c) = -^X^X* with c = n/N > 0. 

Figure 4 plots the density in (32) for different values of c. From (32) the reader may notice that as c ^ 
i.e. for a fixed n as iV — > oo, both the largest and the smallest eigenvalue, and 6_ respectively, approach 
1. This validates our intuition about the behavior of W{c) for a fixed n as TV — > oo. Additionally, from 
Figure 4, the reader may notice that the density gets increasingly more symmetrical about a; = 1 as the value 
of c decreases. For c = 0.01, the density is almost perfectly symmetrical about x = 1. The reader may note 
that with an appropriate scaling and translation, the density of {H^(c)} as c — ^ could be made to resemble 
the semi-circular distribution. In fact, in [14] Jonsson used this observation and the correspondence between 
the distribution in (32) and the generalized /3-distribution to infer the moments of Wn{c) from the even 
moments of the Wigner matrix which are incidentally the Catalan numbers denoted by Ck for an integer k. 
More recently, Dumitriu recognized [15] that these moments could be written in terms of the (fc, r) Narayana 
numbers [16] defined as 

so that the individual moments may be obtained from the moment generating function 

Mf(c) = gc'iV,,. = £c^n('=;l) (34) 

r=0 r=0 \ / \ / 

for which, it may be noted that {!) ~ Ck = M^k are also the even moments of the standard Wigner 
matrix. 



^For notational convenience, when discussing infinite (generalized) Wishart matrices we will simply refer to {W(c)} as the 
Wishart matrix even when the elements of Xn are i.i.d. but not Gaussian. We will occasionally add a subscript such as Wt(c) to 
differentiate between different realizations of the ensemble {VF(c)}. When discussing finite Wishart matrices, we will implicitly 
assume that the elements of Xn are i.i.d. Gaussian. Its use in either manner will be obvious from the context. 
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Figure 5: The density of C„ with c = n/N = 0.1 and di7(r) = 0.6 (S(r - 1) + 0.4 (5(r - 3). 

The Wishart matrix {Vl^(c)} has been studied exhaustively by statisticians and engineers. For infinite 
Wishart matrices, (32) captures the limiting density and the extreme eigenvalues i.e. the region of support. 
The limiting moments are given by (34). As it was for the asymptotic Wigner matrix, the limiting density of 
the asymptotic Wishart matrix did not depend on whether the elements of X„ were real or complex. Similarly, 
as it was for the Gaussian ensembles, the analytical behavior of the finite Wishart matrices did indeed depend 
on whether the elements were real or complex. Nonetheless, the limiting moment and eigenvalue behavior 
could be inferred from the behavior of the finite (real or complex) Wishart matrix counterpart. The reader 
is directed towards some of the representative literature and the references therein on the distribution of 
the smallest [17], largest [18, 19, 20, 21], sorted [19, 22, 23], unsorted eigenvalues [19, 22, 23], and condition 
numbers [22, 23] of the Wishart matrix that invoke this link between the finite and infinite Wishart matrix 
ensembles. We will now discuss sample covariance matrices for which the behavior of the limiting density 
can be best, if not solely, analyzed analytically using the Marcenko-Pastur theorem. 

3.3 Additional examples of sample covariance matrices 

Suppose dH{T) = p6(t — Ai) + (1 — p) (5(r — A2) i.e. r„ has an atomic mass of weight p at Ai and another 
atomic mass of weight (1 — p) at A2. The canonical equation in (29) becomes 



Ai(l — c — cmz) — z A2(l — c — cmz) — z 
which upon rearranging yields the Stieltjes polynomial 

Ai c^TO^z^A2 + (-2A1A2CZ + Aicz^ + 2Aic^A2Z + A2cz^) 

+ (A1A2 + A2C2; + PX2CZ — Xiz + Aic^A2 + — A22; + 2A1C2; — pXicz — 2A1A2 c) m 

— (pA2 + z — pXic + Aic - Ai + PA2C + pAi) = 0. (36) 

It can be readily verified that if p — 1 and Ai = 1, then C„ = -^Tn^^ XnX*Tn^'^ is simply the Wishart matrix 
we discussed above. Though it might not seem so from a cursory look, for p — I and Ai = 1, (36) can be 
shown after some elementary factorization to simplify to (31). Since (36) is a third degree polynomial in m 
it can conceivably be solved analytically using Cardano's formula. 

For general c, p, Ai and A2 this is cumbersome and cannot be solved analytically as a function of z and c 
for arbitrary values of p, Ai, and A2. However, for specific values of c. p, Ai, and A2 we can numerically solve 



the resulting Stieltjes polynomial in (36). For example, when p = 0.6, Ai = 1 and A2 = 3, (36) simplifies to 

3 c^m^z^ + (4 cz^ - 6 cz + 6 c^z) m^ + (3-4z-6c+^cz + 3c^ + z^)m + ^c + z- ii=0. (37) 

5 5 5 

For c = 0.1, (37) becomes 

3 ,0/27 2\ 2 /'243 169 .\ 99 

To determine the density from (38) we need to determine the roots of the polynomial in m and use the 
inversion formula in (5). Since we do not know the region of support for this density we would conjecture 
such a region and basically solve the polynomial above for every value of z. Using numerical tools such as 
the roots command in MATLAB this is not very difficult. 

Figure 5 shows the excellent agreement between the theoretical density (solid line) obtained from numer- 
ically solving (38) and the histogram of the eigenvalues of 1000 realizations of the matrix C„ with n = 100 
and N = n/c = 1000. Figure 6 shows the behavior of the density for a range of values of c. This figure 
captures our intuition that as c ^ 0, the eigenvalues of the sample covariance matrix C„ will be increasingly 
localized about Ai = 1 and A2 = 3. By contrast, capturing this very same analytic behavior using finite 
RMT is not as straightforward. 
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Figure 6: The density of C„ with rfi7(r) = 0.6 (5(r - 1) + QA5{t - 3) for different values of c 

Unlike the Wishart matrix, the distribution function (or level density) of finite dimensional covariance 
matrices, such as the C„ we considered in this example, can only be expressed in terms of zonal or other 
multivariate orthogonal polynomials that appear frequently in texts such as [19]. Though these polynomials 
have been studied extensively by multivariate statisticians [24, 25, 26, 27, 28, 29, 30, 31] and more recently by 
combinatorists [32, 33, 34, 35, 36, 37] the prevailing consensus is that they are unwieldy and not particularly 
intuitive to work with. This is partly because of their definition as infinite series expansions which makes 
their numerical evaluation a non-trivial task when dealing with matrices of moderate dimensions. More 
importantly, from an engineering point of view, the Stieltjes transform based approaches allows us to generate 
plots of the form in Figure 6 with the implicit assumption that the matrix in question is infinite and yet, 
predict the behavior of the eigenvalues for the practical finite matrix counterpart with remarkable accuracy, 
as Figure 5 corroborates. This is the primary motivation for this course's focus on developing infinite random 
matrix theory. 

In the lectures that follow wc will discuss other techniques that allow us to characterize a very broad 
class of infinite random matrices that cannot be characterized using finite RMT. We will often be intrigued 
by and speculate on the link between these infinite matrix ensembles and their finite matrix counterparts. 
We encourage you to ask us questions on this or to explore them further. 
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Exericses 



1. Verified that Silverstein's sample covariance matrix theorem can be inferred directly from the Marcenko- 
Pastur theorem. Note: You will have to do some substitution tricks to get the parameter "c" to refer 
to the same quantity. 

2. Derive the moments of the Wishart matrix from the observation that as c ^ 0, the density becomes 
"approximately" semi-circular. Hint: you will have to make an approximation for the region of support 
while remembering that for any a < 1, < a. There will also be a shifting and rescaling in this problem 
to get the terms to match up correctly. Recall that the moments of the Wishart matrix are expressed 
in terms of the Narayana polynomials in (34) . 

3. Come up with numerical code to compute the theoretical density when dHir) has three atomic masses 
(e.g. dH{T) = 0.4(5(t- 1) + 0.4(5(t- 3) +0.2(5(t- 7)). Plot the limiting density as a function of x for 
a range of values of c. 

4. Do the same when there are four atomic masses in di?(r) (e.g. dif (r) = 0.3 S{t — 1) + 0.25 S{t — 3) + 
0.25 (5(r — 7) + 0.25 5(r — 10)) . Verify that the solution obtained matches up with the simulations. 
Hints: Do all the roots match up? 

5. What happens if there are atomic masses of negative weight in dH{T) (e.g. dH{T) ~ 0.5 (5(r + 1) + 
0.5(5(r— 1)). Does the limiting theoretical density line up with the experimental results? Check the 
assumptions of the Marcenko-Pastur theorem! Is this "allowed" ? 
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